
#include "problem.h"

class Problem1 : public Problem
{
public:
  Problem1() ;
  
  ~Problem1() {} ;
  
  // boundary conditions
  bool DirichletBC(const Point& pt, Real& val) const;
  
  bool NeumannBC(const Point& pt, Point& val) const;
  
  bool RobinBC(const Point& pt, Real& val1, Point& val2) const;
  
};

Problem1::Problem1() : Problem(2, 1, 1)
{
  const double PI = 3.141592653589793 ;
    //GetPot input_file("/emchome/vijaysm/KARMA/input/2d_homogeneous_eigen_problem/neutronics_simple_eigen_input.in");
	
    //const double len = input_file("length", 1);
    //dimension    = input_file("dim", 1);
  const double len = 10.0 ;
	
	lx = ly = lz = len ;
	
  // 1 group cross-section data
	xsrem[0][0] = 1.099822E-02 ;
	xsfiss[0][0] = 0.4106224E-02 ;
	xsnufiss[0][0] = 0.0110556 ;
    //nu = xsnufiss[0][0] / xsfiss[0][0] ;
	xsdiff[0][0] = 1.360958/100 ;
	xsscatt[0][0][0] = 0.0 ;
	chip[0] = 1.0 ;
	
	double Bg2 = PI*PI/lx/lx + (dimension>1)*(PI*PI/ly/ly) + (dimension>2)*(PI*PI/lz/lz) ;
	
    //double nu = 2.6924006094163396833684669905977 ;
	
	kinf_exact = xsnufiss[0][0] / (xsrem[0][0]+xsdiff[0][0]*PI*PI/lx/lx) ;
	
	PetscPrintf(libMesh::COMM_WORLD,"\n%D-Dimension independent 1-group eigen problem with Keff =  %G\n", dimension, kinf_exact) ;
}


#undef __FUNCT__
#define __FUNCT__ "DirichletBC"
bool Problem1::DirichletBC(const Point& pt, Real& val) const
{
	const Real x = pt(0), y = pt(1), z = pt(2) ;
	
	libmesh_assert(dim <= 3) ;
  
  bool success = false ;

	if (dimension == 1)
		if(fabs(x) < 1E-12 || fabs(x) > lx-1E-12 )
		{
			val = 0.0 ;	// Or do nothing, since by default, initialized to 0.0
      success = true ;
		}

	if (dimension == 2)
		if((fabs(x) < 1E-12 || fabs(x) > lx-1E-12 ) || (fabs(y) < 1E-12 || fabs(y) > ly-1E-12 ))
		{
			val = 0.0 ;	// Or do nothing, since by default, initialized to 0.0
      success = true ;
		}

	if (dimension == 3)
		if((fabs(x) < 1E-12 || fabs(x) > lx-1E-12 ) || (fabs(y) < 1E-12 || fabs(y) > ly-1E-12 )
				|| (fabs(z) < 1E-12 || fabs(z) > lz-1E-12 ) )
		{
			val = 0.0 ;	// Or do nothing, since by default, initialized to 0.0
      success = true ;
		}

	return success ;
}

#undef __FUNCT__
#define __FUNCT__ "NeumannBC"
bool Problem1::NeumannBC(const Point& , Point& ) const
{
	return false ;
}

#undef __FUNCT__
#define __FUNCT__ "RobinBC"
bool Problem1::RobinBC(const Point& , Real& , Point& ) const
{
	return false ;
}

